
# 
# This script reproduces and applies code from the replication files of Yimeng Li (2019)
# "Relaxing the No Liars Assumption in List Experiment Analyses"
# DOI: https://doi.org/10.1017/pan.2019.7

remove(list=ls())

# Load libraries:
library(bindata) # required by sim_bounds.R
library(dplyr)
library(forcats)
library(haven)
library(ggplot2)
library(list)
library(here)

library(nleqslv) # required by list_relaxed_liars.R
library(purrr)
library(tidyr)

library(patchwork)
source("Li_2019_NoLiars/list_relaxed_liars.R")
source("Li_2019_NoLiars/list_types_plot_US.R")
source("Li_2019_NoLiars/list_create_df.R")

list_dataUS <- read_dta("data/listUS.dta")
list_dataUS <-as.data.frame(list_dataUS)
list_dataUS<- list_dataUS  %>% 
  mutate(treatfac=as.factor(treated),
         IDvar=as.factor(IDvar),
         round.fac=as.factor(round)) |> 
  mutate(value= value)

options(digits = 3)


dist_control_Example <- c(151, 366, 909, 107, 38)
dist_treated_X_Example <- c(91, 247, 527, 592, 74, 41)

DiM_X_Example <- ictreg(y ~ 1, data = create_list_df(dist_control_Example, dist_treated_X_Example, J = 4), treat = "treat", J = 4, method = "lm")


unname(DiM_X_Example$par.treat)
unname(c(DiM_X_Example$par.treat+qnorm(0.025)*DiM_X_Example$se.treat, DiM_X_Example$par.treat+qnorm(0.975)*DiM_X_Example$se.treat))

bounds_X_Example <- list_relaxed_liars(dist_control_Example, dist_treated_X_Example, J = 4, affirmative = TRUE)

bounds_X_Example$bounds
bounds_X_Example$CI

FigureUS<- list_types_US_plot(bounds_X_Example) 

F1<- FigureUS$plot_two_types
F2<- FigureUS$plot_three_types

F1+F2
ggsave("FigureA11.png",
       path = here("figures"),
       dpi=600)
